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ABSTRACT 



With the advent of the high speed digital computer, many 
problems heretofore considered unsolvable for all practical 
purposes are now well within the reach of the aoplied math- 
ematician. One such problem is the routing of a ship through 
a time dependent ocean wave field, from one point on the 
earth's surface to another, so as to minimize a cost func- 
tion of the form g(x,y,t,u). 

This paper considers a numerical solution to the above 
problem. The technique to be employed is known as the method 
of steepest ascent and is attributed to Arthur E. Bryson and 
Walter F. Denham [1]. Although the computer program as given 
in the Appendix is written specifically for a VC2AP3 class 
vessel operating in a described area of the North Pacific 
Ocean, it can be readily modified to accommodate any type 
vessel operating in the Northern Hemisohere. 
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I. statement of the Problem 

The basic problem to be developed and solved in this 
paper is the following: A VC2AP3 vessel, located at latitude 

40°N longitude 154°E, is to transit a time dependent ocean 
wave height and wave direction field to latitude 38®N longi- 
tude 123°W so as to minimize a cost function of the form 
g (x,y ,u,t) . 

A south-polar stereographic projection of the Northern 
Hemisphere upon a plane passing through the circle of 60°N 
latitude is used to establish a rectangular coordinate sys- 
tem, with the OX and OY axes parallel to the projections of 
the meridians of 10°E and 100°E longitude, respectively. 

Then a 63 x 63 grid is constructed along the OX and OY axes, 
with X = y = 31 defining the projection of the North Pole. 

At each point of this grid, actual wave heights and wave 
directions taken from the files of the U.S. Navy Fleet Numer- 
ical Weather Central, Monterey, California, are recorded. 

The map scale factor MSF, defined as the ratio of a differen- 
tial distance in the OXY plane to the corresponding differen- 
tial distance on the earth's surface, is 

MSF = [973.75+(x-31)^+(y-31)^]/1043.6. 

Germane to the solution of the above problem is the 
elliptical polar velocity diagram of Figure 1, giving the 
ship's speed v as a function of the angle 0 between the 
ship's heading and the wave direction. R. W. James [6] gives 
empirical curves for the speed of a VC2AP3 class vessel in 
head waves v^^ , beam waves Vj^ , and following waves v^ as a 
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function of wave height. Bleick [2] showed through a least 
squares analysis that these three speeds can be represented 
closely by the 4-parameter formula 



V = C^(H+C^) ^ exp(-C3H) , 

where H is the wave height. 

Thus, the projected speed of the ship at any point on 
the grid is 

V(x,y,0,t) = v(H,0)MSF , 

where v(H,0) is the ship's speed of the polar velocity dia- 
gram of Figure 1, and H(x,y,t) is the interpolated value of 
the Fleet Numerical Weather Central wave height data. If 
we define K(x,y,t) to be the stereographic grid wave direc- 
tion measured counterclockwise from the OX axis, and u(t) 
to be the ship's heading as measured counterclockwise from 
the OX axis, then 

[v sin(u-K)/b]^ + { [c+v cos (u-K) ]/a}^ =1 or 
v(H,0) = ^/[— COS0 + /b^cos^0+ (a^-c^) sin^0] 

3i ^ 

where 0 = u-K, a and b are the semi-principal axes of Figure 
1, and c is the distance to the eccentric pole o. We note 
at this point that an elliptical polar velocity diagram, 
such as Figure 1, must be specified for each wave height, 
and that a, b, c are functions of H(x,y,t). Thus, we can 
write the equations of motion of the ship's projection in 
the OXY plane as : 

x(t) = Vcos u(t) 
y(t) = Vsin u(t) 
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Figure 1. Polar Velocity Diagram 



where the ( 
and V = |v| 
the control 
to minimize 



) over a variable indicates its time derivative, 
= V(x,y,u,t). Hence, our problem is to choose 
angle, u(t) , at each point of the transit so as 
the cost function g. 
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II. Abstract of the Theory 

Consider a system of differential equations of the form 

= f^ (i=l,...n), where f^ = f^(X,U,t) are functions of 
2 

class C , X is an n-dimensional vector, U is an m-dimensional 
vector, and t, the independent variable, is a scalar. We 
desire the solution to the above system of equations which 
is optimal in the sense of maximizing (minimizing) the ter- 
minal value of some cost function g(X,t)^_^, while simulta- 
neously satisfying terminal constraints, 

= 0 (i=l,. . .N; 0 < N < n) . 

The solution of a problem of the above form by the method of 
steepest ascent is concerned with the development of the m 
control variables, u^^ ( t) , . . . u^ (t) . Essentially the problem 
is solved in two parts. In part I we are interested in 
attaining admissibility, that is satisfying the terminal 
constraints 

h^(X,t) = 0 (i=l,. . .N; 0 < N < n) . 

In part II we are concerned with the problem of maximizing 
(minimizing) the terminal value of the cost function g. The 
procedure is as follows: We initially guess nominal values 
of the control, U; in general, the resulting solution will 
neither be admissible nor an extremal. Changes in the con- 
trol are then generated which tend to drive us to admissibil- 
ity, that is, U(t) is replaced by U(t) + 6U(t). This nro- 
cedure is continued until we attain admissibility. The ad- 
missible solution, in general, will not be an extremal; so 
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we once again generate changes in the control which tend to 
drive us toward an extremal without affecting admissibility. 

The author realizes the above is rather abbreviated; 
however, this discussion is intended to give a brief idea of 
what is to follow. 
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III. Development of the Theory 

Notation 



Let us start by considering the following system of 

ordinary nonlinear differential equations: 

•1 • „ ^ 1,12 
X = X = V cos U = f (X ,X /U,t) 

(3.1) x^ = y = V sin u = f^ (x^ ,x^ ,u,t) 

= z = g(x,y,t) = f^ (x^ ,x^ ,u, t) 
where V = V(x,y,u,t) , u = u(t) , and g is the cost function 
to be minimized. We can write this as 



x^ = f^(x^(t) , x^(t) , u(t) , t) , (i = 1,2,3) , 0 < t <T, 
where t, the independent variable, is a scalar. It will be 
convenient in the development to follow to represent these 
sets of variables by matrices or by vectors , so let us 
define : 





(t) 




f^ (x^ ,X^ ,U,t) 


->■ 

X = 


x^(t) 


II 


f^ (x^ ,X^ ,U, t) 




x^ ( t) 




f^ (x^,x^ ,u,t) 



where U is a single-rowed column matrix since there is but 
one control variable. 

Equations (3.1) then become 
(3.2) X = F 

The vector X is called the state vector , or state , and the 
one dimensional vector U, the control vector , or control . 
We will assume the control to be a continuous function of 
time . 
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Now the problem as posed in section I can be intemreted 
as follows: We desire to find a curve C, a solution to (3.2) 

for some choice U, which is optimal in the sense of minimizing 
the terminal value of the cost function g while simultaneouslv 
satisfying the terminal constraints 



x^(T) = x^, x^(T) = xj , 



1 2 

where x^ and x^ are the final values of the state corresoond- 
ing to the rectangular coordinates of the terminal point in 
the OXY plane . 

The Adjoint System of Equations 
Of primary interest is the relationship between varia- 
tions, 6u(t) , in the control at any time, t (0 _< t ^ T) , and 
the resulting terminal variations of the state variable, 
6X(T). If we consider (3.2) in component form, the varia- 
tional relationships for the state variables can be expressed 
as : 



(3.3) 



^ -1 
ox = 



x*2 

6x = 



3f^ 


X 1 ^ 


1 — 1 
m 

CO 


^ 2 , 


3f^ 


r 3 , 


1 — 1 
M-l 
CO 


6u 


a 1 
3x 


6x + 


3x2 


6x + 


3x2 


6x + 


3u 


3f^ 

1 


^ 1 ^ 
6x + 


3f2 

0 


6x2 


3f2 


X 3 ^ 
6x + 


3f2 


5u 



9x 



9x 



3x' 



3u 



.•3 3f^ . 1 ^ 9f^ . 2 ^ 9f^ 3 ^ 3f^ ^ 

ox = :p OX + j OX + o OX + ou 

3x"^ 3x'^ 3x 3u 



Writing this in matrix form we have: 



6X = G6X + H6U 
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where 


9f^ 


3fl 


3f^ 




'afi' 


1 






dx^ 


3x2 


3x2 




9u 




• 1 
















6x 


G = 


dx^ 


3x2 


3f2 

3 

3x 


1; H = 


3f2 

3u 


• 

; 6X = 


6x2 

6x2 






3f2 


3f2 




111' 










3x2 


3x2 




3u ' 


1 

1 





and 6U = [6u] . 

— y* 

Note in the special case where the control U consists of a 
single variable, u, 6U can be considered a single-rowed 
column vector. 

If we now introduce a new set of variables 

X(t) 



A(t) = 



VI (t) 
a(t) 



Lagrange multipliers, which are undefined functions of t, and 

T T 

multiply (3.4) through by A , where A is the transpose of 



A, we obtain the scalar relationship 

(3.5) a'^6X = a'^G6X + a'^H6U . 

Let us now integrate (3.5) over the interval (0,T) to obtain 



(3.6) /q a"^ 6X dt = /J a"^G 6X dt + a'^H 6U dt . 



Integrating the integral on the left by parts yields 

a"^6X|q - /q a"^ 6X dt 
Rewriting this as 



m m ^ m m ^ 

A G6X dt + /J: A H6U dt 



a'^ 6x|q - /q (a”^ + a'^G) 6X dt = /q a'^H 6U dt 
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and choosing A as a solution to (A*^ + a'^G) 6X = 0 defines the 
system of equations adjoint to the variational equations (3.4). 
That is, the adjoint system of equations is 



(3.7) 

With this choice of A 

(3.8) = 



A^ = -A^G 



(3.6) reduces to 



T I rT T 

A 5X| Q + /q a H6U dt 



and, since we are assuming 5X(0) = 0 , we have a relation be- 
tween the terminal variations of the state variable, 6X(T), 
and the variations in the control, 6U. Now if A is chosen 
as the specific solution to the adjoint system of equations, 
(3.7), such that at time t = T, A(T) = (1,0,0), then (3.8) 



becomes 



6x^(T) = /q a'^H 6U dt 

giving the terminal variation of x due to variations 6U(t). 

th 

Similarly, choosing A^(T) such that the i comoonent is 1 
and all other components are 0 yields 

(3.9) 5x^(T) = /q aTh 6U dt 

T 

We may consider the matrix product A^H as a vector, say A^^(t), 
when there is more than one control variable. Then (3.9) can 
be written as 



5x^(T) = /J A^(t) • 6U(t) dt 

and for a maximum change in x^(T) we would choose 6U(t) par- 
allel to A^(t) , for 0 £ t £ T. Courant [3] oaqe 223, refers 

“> i i 

to the vector A^(t) as the gradient or x , of A^ = Vx . 
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Special Variations 



From the previous discussion we have seen for a maximum 

i ^ i 

change in x we would like 6U(t) to be parallel to Vx for 

0 _< t _< T. This motivates the following choice of a soecial 

variation of the control, 

(3.10) 6U(t) = ejL^x^(t) + e2^x^(t) . 

It might be worth noting we do not consider ^x^ since at this 
juncture we are only interested in satisfying the terminal 
constraints 

l/^^ 1 2,.. 2 

X (t) = x^, X (t) = x^ . 

T i 

Substituting (3.10) into (3.9), and recalling A^H = Vx , we 
obtain 



X 1 

fix 


(T) = e 


J 0 




• Vx^ 


dt + e2 


rT 
J 0 






dt 


11) 


^ 2 
fix 


(t) = e^^ 


j 0 




• Vx^ 


dt + e2 


j 0 


^x2 


• ^x^ 


dt 



We now define 

= /q Vx^ • Vx^ dt 

so that (3.11) becomes 

2 

6x^(T) = I e.Z^^ (i=l,2) 

j = l ^ 

For small perturbations, the value of T will be changed 
only a small amount dT, so that 

dx^(T) = 6x^(T) + x^(T) dT (i=l,2) . 
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We can express this relationship in matrix form as 



r 1 

1 ^11 ^12 •! 




- 


i 


1 

1 

I — 1 

1 


I Z Z X 




®li 


i 

1 _ 


X^ - X (T) 


' „21 „22 *2 








2 2 


Z Z X 

L. 


1 


®2 




x^ - X (T) 



, dT 

or ZE = AX. Let us now digress for a moment, and recall that 

if we have an overdetermined system of ecruations , say AX = B, 

where A, the coefficient matrix, is 3x2, X, the unknown 

matrix, is 2 x 1, and B is a 3 x 1 matrix, then a least 

squares solution for X is obtained by multiplying both sides 

T T 

of the above equation on the left by A provided AA is not 
singular. Since we have a system of two ecruations in three 
unknowns, let us apply the above technicrue in reverse. 

Let E = z"^C where C = 

T T 

becomes ZZ C = AX. Now ZZ is a 2 x 2 Gram matrix, so let 

T 

ZZ = A = 

Thus (3.12) is representable as AC = AX, and we are now able 

T 

to solve for C. Use of (3.10) and E = Z C gives us our var- 
iations in control 6U(t) , for 0 £ t £ T. 

At this point in the development we should be able to 
obtain an admissible curve. We will now consider the prob- 
lem of decreasing the cost function V'?ithout affecting x^(T) 

2 

and X (T) . To accomplish this end we consider a special 
variation of the form 



11 


^12 


21 


^22 



so that (3.12) 
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(3.13) 6U(t) 



= aj^^x^(t) + a2Vx^(t) + (t) 



Substituting (3.13) into (3.9) yields 



6x^(T) = ^ a.Z^^ and 

j=l ^ 

dx^(T) = 5x^(T) + x^(T) dT 
In matrix notation this becomes 



(i=l,2,3) 





1 

^f 


1 

X 




z^i 


Z"2 


zl3 


x^(T) 


(3.14) 


2 

^f 


- x^(T) 


= 


z21 


z22 


z23 


. 2 

x^(T) 




. 3 


- x^(T) 




z^i 


Z^2 


z33 


x^(T) 


Since 


we 


have an 


admissible path. 








1 

^f 


1 

- X 


(T) 


2 

^f - 


x^(T) 


^ 0 



a. 



a. 



a 



3 

dT 



3 3 

however, x^ - x (T) is un)cnown, and, in general, is not zero 
Once again we are faced with the problem of solving an under 
determined system of equations, (3.14). As before we will 



let a = Z D where D = 









so that (3.14 becomes AX = ZZ D. 



T 

Now ZZ is a 3 X 3 symmetric Gram matrix so let 



(3.15) 





^1 


*^12 


*^13 


T 

ZZ = B = 


*^12 


*^22 


^23 




*^13 


*^23 


^33 


is representable 


as 






AX = 


BD 
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Define = £3^^ (i=l/2,3), where e is to be determined, 

and cofactor of in the determinant of B, that 



IS 



di = e [b ^2 ^=*23 


^22 


CO 
1 — 1 


d 2 = e [b ^2 *=*13 


1 — 1 
1 — 1 

1 


^23^ 


d 3 = e [b^^ b 22 


^12 


^ 12 ^ 


Then, upon substitution into 


(3.15) 


/ we 







X j - x^ (T) 




0 


(3.16) 


AX = 


Xj - x^(T) 


= e 


0 






x^ - x^(T) 




det B 



Thus, by choosing e negative, we can decrease x and keep x 
2 

and X fixed so long as B is nonsingular. That is 

12 3 

dx dx „ , dx'’ ^ 4 . T. 

• 3 — = - 3 — = 0 and 3 — = det B. 
de de de 

T -*■ 

Now, having D and a = Z D, (3.13) may be solved for 6 U as a 

multiple of e. In order to limit the maximum change in the 

control, we define the following norm: 

dT[ 



N = {max u (t) + 1 



} 



We now define e = - 



t 

D ELUMX 
N 



, where DELUMX is the maximum 



allowable change in the control. Thus 6 U(t), the solution 
to (3.13) is replaced by 

e 0 U(t) 6 U(t) 
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IV. The Numerical Solution 



The numerical solution of the problem as posed in Sec- 
tion I is based upon the solution of the original system of 
equations, the adjoint equations, and the variational equa- 
tions. For the sake of simplicity, a finite difference 
approximation to the derivative with a variable time step of 
integration is used to integrate the above equations. The 
variable time step of integration is obtained by dividing 
the total time as determined after each path into two hun- 
dred equal increments. Initially, the value for T, the 
terminal time, is determined by the geodesic track, but this 
is not a requirement. Any reasonable guess for T will suf- 
fice. Justification for this method of integration is based 
upon a comparison of the geodesic route as obtained, using 
the above mentioned technique and a fourth order Runge-Kutta 
integration method. 

Due to the core storage limitation imposed by the Fleet 
Numerical Weather Central, 165,000 octal words, a theatre of 
operation is established which is a subfield of the 63 x 63 
grid mentioned in Section I. Having established this sub- 
field the wave field data is then read into core memory. We 
are now in a position to proceed with the numerical solution 
of the problem. 

Although not a requirement, the geodesic track from the 
point of departure to the point of destination is calculated 
for comparison purposes. As stated in Section II, a nominal 



19 



control program is initially required. This is taken as the 
angle of inclination of the straight line connecting the 
track end points. 

Recall that in the development of the theory we chose 

as particular solutions to the adjoint equation, A 2 , and 

th 

A^/ where A^^ (T) has its i^ component equal to 1 and all other 
components equal to 0. To satisfy this requirement we define 





1 

1 — 1 

L 


1 


0 




1 

0 

1 


> 

0 

II 


0 0 

1 


, A2(0) = 


1 

0 


> 

U) 

0 

II 


1 

0 fH 



and integrate the adjoint system of equations forward to 
obtain A^^ (T) , A 2 (T) , A^(T) . The resulting matrix 



Xj_(T) 




Xjd) 




P2<T) 


M3(T) 


1 — 1 

D 


02(T) 


03(1) 



is then inverted to obtain new initial values A^(0), A2(0), 
A^CO) so that upon integration forward a second time one 
obtains , ' « 



A^(T) = 0 , A 2 (T) = I 1 I , and A 3 (T) = 



0 
I 0 
1 ! 



Having obtained the initial values of the adjoint vectors 
Aj^, A 2 , and A^ such that at the terminal time T they take on 
prescribed values, we will run a oath, check for admissibility, 
and if not admissible generate corrections to the control. 
Simultaneous integration of the original system of equations 
and the adjoint equations gives a path. At each point along 
this path the gradients in the direction (i=l,2,3) are 
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calculated and stored for later recall in determining the 
changes in the control. In addition, the elements of the Z 
matrix of Section III are calculated. A check is now made 
to see if we have attained admissibility. If not, we pro- 
ceed to calculate the changes in control which will tend to 
satisfy the admissibility constraint. Assuming we do not 
have an admissible path, the A, C, and E matrices of Section 
III are calculated. Thus we obtain new values for the con- 
trol and the terminal time T, that is, U(t) is replaced by 
U(t) + i5U(t) and T is replaced by T + dT. The above proce- 
dure is continued until the admissibility constraint is met. 

Having admissibility we consider the problem of de- 
creasing the cost function 

fT Y H,. 

g = /q e’ dt , 

where y is a scalar, without affecting admissibility. A 
check is made to see if the terminal value of the cost func- 
tion as determined by the last admissible path is less than 
that determined by the previous admissible path. If it is, 
we proceed to calculate the changes in control as defined in 
Section III, so as to decrease the cost function g. If not, 
we decrease the maximum allowable change in the control at 
any time t by one-half. A check is then made to determine 
if the maximum allowable change in control is less than 0.01. 
If it is, a solution has been obtained. If not, proceed to 
calculate the changes in control so as to decrease the cost 
function. Assuming we have no solution we calculate the B, 

D, and a matrices of Section III. This enables one to obtain 
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new values of the control, that is, U(t) is replaced by 
U(t) + 6U(t). A new path is run which, in general, is not 
admissible; so we again drive to admissibility. This oroce- 
dure is continued until a solution is obtained. 

It should be noted the values used in the program for 
the maximum allowable change in the control and the stopoing 
criterion are quite arbitrary. 
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V. 



Conclusions 



The value of this study lies not in the solution of the 
particular problem considered, but in the fact that given 
any realistic cost function which is a function of the 
stereographic grid coordinates x and y, the control U, and 
time, one can obtain an optimum solution with only slight 
modifications. Such a cost function for a cargo type vessel 
might reflect costs due to storm damage, failure to meet 
scheduled commitment dates, and spoilage of cargo due to ex- 
tended time at sea, to mention but a few. 

On the basis of the results obtained the following 
observations seem justified: 

(a) It was shown in Section III that a neces- 
sary condition for the solution curve C to be an ex- 
tremal is that the determinant of the B matrix at 

t = T approach zero as track iteration continues. 

In the .case where we chose y = 0, that is, minimum 
time, we saw that this, in fact, was the case. In 
the second problem, where we chose y = 0.01, the 
determinant of the B matrix was reduced from a value 
on the order of 5,000 to 0.9. The conclusion to be 
drawn here is that, although the solution curve C 
is not an extremal, it very closely approximates one. 

(b) Defining the initial control as the con- 
stant function u(t) = V (0 £ t £ T) , where V is the 
initial angle of departure, leads to convergence to 
within fifteen miles of the point of destination. 
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Thus there seems to be no problem as to how to 
initially choose the control so as to obtain con- 
vergence . 



(c) The simplified integration routine used 
is feasible and yields satisfactory results. 

(d) The requirement of storing an entire con- 
trol program and a complete set of functions 

T 

i = 1,2,3 for any one iteration is a marked 
disadvantage. Due to the large storage reauire- 
ments imposed by the establishment of the wave 
fields a relatively small amount of core memory 
remains . 

(e) An analysis of the results obtained show 
a substantial decrease in the wave heights for the 
minimum time route, as compared to the geodesic 
track, with a decrease in the total transit time 
of approximately 2% hours. The weighted time 
route, on the other hand, shows only slight de- 
creases in the wave heights encountered, as com- 
pared to the minimal time track. This leads one 
to conclude that if the objective is to limit the 
wave heights encountered, a function other than 
exp yH should be used. This is more readily seen 



if we expand e'' , that is, 

e^^ = 1 + vH + + . . 
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For small values of y and nominal wave heights, H, 

/*T yH 

the dominant factor in minimizing z = Jq dt is 



time . 
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APPENDIX 



COMPUTER PROGRAM 
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( IMPiir»OUT^’nT»TaPE^) 

' r .' t NSi OfM "Tr; ( 1 ? ) ,ii (?no ) ,nF:i.ii (200 ) »gpadx (200) ,nRADY (?oo ) « 

1 rPS(? 01 )»HT( 201 )tAZ{ 2 nn»AL(? 0 l)«WA( 20 l) ♦ wE ( 39 H 9 ) »Wn ( 63 * 63 ) 

F. v|.AM( 3) *TI AM(3) *XMU(3) tXSK-iMAO) *fiPA07(20 0) *TMII(2) 

COMMON! v,Y»A*n,CC*H*CK,SK,KX*KY*KT,FOXt FO Y * KCON 

COMMON/l l/XhT(044P) ♦ SNK (044 0) *T 

fOMMONj/L 2/HX ,HY 

C'1MmOM/L0/CAPV ,CAP\/X , C APV Y ♦ CAPVU » ia> 

^ r OMM0N/L6/XK . XLG* XLT 

t (.Mil VALFmCE {-'F ( 1 ) *w1M I * 1) ) 

►vF AL M0P 

I l*'Tt(^FP '■' 1(3 

^ /'f 1- II L FIPST ESIARI TSH THF '*'AVF FTEur) IM OUo THFATFP OF INTFOPST 

I TrAUXsn PT0X = n ^ IMIMXsO 'S J.MIMX = 0 'f> KTX = 0 

^ pf'M) otmfNSTON'; AMn 'ijnijmiim TNJnICp'5 Of apray0 
p(-' A U 0 , H V , K Y . K r * r M I M ♦ JM T M , CR T T 
■n ^np^lAT ( 0 T 2 .F 4 .O ) 

I o-’Iimt 3.KX ,KY *kT ♦ jmtn ,r.RlT 

I; Y F opmaT ( 3'^hMOTMi: M«;T 0 M Oh WAVE field APR A YSr (1 3 * 1 H * T 3 ♦ 1 H ♦ 1 3 , 1 H ) 4 X5 w T 
lMlM=r ^,4y0HjMlN=r3.4x5HCRlT = F4.O ) 

,r>v = 32-J*i1 i-.i 

! FOy=I()x 

|: ( OY=JOY 

'' h c 0 M = K X K Y ♦ K X 

x^-F AMsFI OAtI- (Kx* 1 ) /?, 
y.iFAcj = M OaTF (KY*1)/2. 

U -»X = XmF AN'“2. a 1 
POysy'^iF 1 

p Cl La*' ('TO aRpay 

,iO b I =1 ,K T 
f' I . T o ( I ) = 0 

' kh 'i<j octal I'ATF.riME OROHO'^ OF THF Tl'.^F PFPlFR MFMRFR0 
kh dij ) , ii T (’ ( I ) * T 1 AU 

I 1 I (•■•■'Pf3'\T(n20«^x, I?) 

I I T ( dTGX .F 0. i)TC, ( 1 ) , aMO. T TaIIX .FO. ITa()« aMI). TMiMx ,F0. IMIm, AND. 

) .iMJiMX ,F 0, J'^'TM. AMji.K r Y .fO.KT ) 750 *620 

SF I array oTG to SE'-FcTEI.I Tiff SFPTES 
i J ' = K T-J 

I ' T 1=1 , JJ 

( .>LL f.'TC^h ( DYG ( T ) . ITAII.DTG ( I *• 1 ) ) 

7 0>M T T Ml iF 

p pt., THE liMHiiCKfD aRRaY«; 

PF..ilNf) 4 
I ) 7 0 n K = 1 , K r 
KKK=( ( K-1 )*KY-l )»KX 

MO HSn KK=) »? 

j Ni 0 = U 

K.*i = KK-1 

c ALL F TFI 0 ‘^X ( W'F .OTG (K ) ,K w,NO) 

IF (^0-l )90,FW,9n 
Hw PPINf Ql .OTOIK) 

01 F OPkAT ( 1 4H()F lELO OF DTO = o20 » 1 X31HNOT FoUMt) ONI INPUT TAPF4 , ARoRT/) 
GO U) 0 3 

Rn GO TO (P 2 » 04 ) KK 

02 IMy H 1=1 ,K X 
1 I=KKK*I 
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()Fi..x = i-Toy 
DO H Jsl.KY 
JJ=KX*J*II 

ofly=j-joy 

DonT = SOf^TF ( OFLX<»nFLX ♦ DEI.Y*OF.LY ) 

( T ♦Tm1m,J*.JMIM) /«;, 729577951 
COSsCOSF ( ARO) 

S» WsFTNF (ARO) 

('9k ( JJ) s (-nELX«>COS-DFL Y*SIN) /root 
SMK (vJJ) = (OFL X«STN-DFLY«COS) /ROOT 
(^n 70 

V4 I'O in 1 = 1 ,KX 
TT= kKK + I 
DO 10 J=1,KY 
J.J = KX*J + IT 

in XH7 (JJ)=wn(I*lMIN,J*JMI(M) 
iS50 COMTTNUF 

IF(SFMFF S/'irCH l)A,lr),700 
Ain POINT Iloo* K,OTO(K) 
lion format () H n» I?, 2 H. ,o?n) 

DO 620 Jsl.KY 

A/() POINT 1?00. ( WH( I + THIN, J4JMTN ),I = 1»KX ) 

1200 fORMAK 1 x»^?6F5,1/ ) 

700 conTTNOF 
PF(«Imo 4 

DTOX = mG(n 'S ITAllX = ITAtJ f. ImtNX=TMTM Ti JMINX=JMTM $ KTXsKT 
OO 750 Tsl.KX 
i'O 750 Jsl.KY 
jJsKXOJ .. I 

RT=I*IMIM % 0J=J+JMIM % KF«=0 
call SFAIAND (RI.RJ.LAND.KER) 

IF (KER) 546,640,646 
64A PRINT 647 S GO TO 53 

64^ FORMAT (26H0RI ,RJ INDEX ERRORS, APORT/) 

A/.« IF (LAND) 649. 750,649 
649 JO 749 K=1,KT 

KKK= ( (K-1 ) «-KY-l ) «KX ♦ JJ 
740 XhTJKKK 1=20.0 
750 C.OM7INUF 

r '«f Should now have dur wave field daTa stored in mfmmory 
<- RFaD hour of DERaWT'JPF and coordinates Of track end points 
RFAO 13.HR,XL61 , XLTI .XL62»XLT2 

1 3 format (F3.0,F6,1,F5.1 ,F6,1 ,F5.1) 

PRINT 14,hR«XLG1»XlT1.XLG2»XlT2 

14 F0PMAI(21H ROUTE OF SHIP BEG I NSF4 , 0 ♦ 3 1 H HOURS AFTER TIME SERIFS OR 
1IG1N/19H fFOM longitude * F6.1»16H AMD LATITUDE = F6.1*/19H TO 
2 I.ONGITHDF = F6.1,16H ANQ LATITUDF = F6 , 1 ) 

r 9F mill now convert coordinates OF TRACK EnD POINTS TD GPlO VALUES 
ARG= ( XLGl-1 0, 1 /5 7, 29677951 
COSLGl=COSE ( ARC) 

SIMLG1=SIMF(ARG) 

ARG=(X| G?- 1 0, ) /57. 29577951 
C0SLG2=C0SF ( ARG) 

STnLG2=STNF ( ARG) 

A PGsXLT 1/5 7. 29577951 
COSLTl=COSF ( ARG) 

SIULTlsSlNF < ARG) 

ARG= XL TP/57. 29577951 
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C0SLT2=C0Sf- ( AWG) 

«;tml T?=«;IMF(AWG) 

.?n'S<H:0'?LTl / ( 1 . ♦SINLT 1 ) 

X 1 =ppl*»rnSLGl 

Y l=PPl*SlNiLGl 

PP? = 11 ,?n5«COSL.T2/ ( 1 .♦SIMLT2) 
n?=Pp?<»rnSLGP 

Y P = p»p*c;] M( G? 
v';T«WTsY 1 ♦Ff’ V 
YGTAPTsYI *F0Y 
vF ^,r)=x2+Fox 

> F M I ' = Y 2 ♦ f f' Y 

r fOtnTF GFOPfSlC PODTF 

M =SINl. T>>«>rOPl.Tl *STMl Gi - COSL I N| T 1 *S T ML G? 

f.M= - STNl.T2<»rOGi.,Tl*COSLGl ♦ C‘'SL T?<>S I MLT 1 ^'’COGl G2 
t'ls ( SIN'l G2*C0SLG1 - C0SLG2*«;INLG1 ) <»CO«;i T1 *cn‘>LT2 
POOTsSG'PTF (FL<»FI. ♦ FM^'FY > EM«FW ) 

M.=EL/R0OT 
F ''rtM/POni 
t ri=EM/P00T 
nf LX = x?-x I 
I Fly-YP-y 1 

Sl? = SQRTF (PFi.XWRFl.X > nF| Y»|)FI Y) 

^PC=S12 

JF (XLGP-Xl-'-] )20,?1 »20 
!\pG = ABSF (FM/f,;?.41 ) 

:\PC = AS1IjF ( aRG<>S1 2) / APG 
21 >=xSrART 

Y = YS TART 
T-HR/P4. 

. TFSrO.n 

f. all Tf RP 
nT ( 1 ) sh 

1 F S ( 1 ) = X T F G 
CALL AMGLF 
A7 ( n SXL.A; 

..I ( 1 ) =Xl T 

«; A ( I ) =X,K 

MK=1 

)f X = XSTA!-'T 
Y YsYSTART 
SG=0.0 
l,.R = 0 

STFP=1 ./F.. 

■)0 40 K=?,?0() 

C.'''SP= (FOY-Y) «FN/31 .?0b + EM 
4TN)P= ( X-FOX ) ^^F^J/31 ,205 - EL 

call TFRR 

CALL ShTR 

CALL VOFRiw 

x = XX ♦ CaPV»COSP«»FTFP 

Y = YY ♦ CApw<>3lNP*SrEP 
S=SS ♦ CAPV*STEP 

IF (S-APC) 35,34,34 
34 P«T= ( APC-SS ) / ( S-SS ) 

t=rat<*stfp *■ T 

XTFS=RAT«STFP + XIFS 
X=(X-XX)4PAT ♦ XX 
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Y= ( Y-YY ) ] ♦ YY 

GO TO 15 
3b XV=X 
YY = Y 

T=T ♦ STFP 
XTFS=X1FS ♦ STFP 

,1F( (K-l)#2/6 ♦ 1 - MK)1(S, 16,15 

]b fuK=MK * 1 

call TtPP 
MT (NK ) =H 
TPS (MK)sXTPS 

call amglf 

A7 (NK ) =yl.G 
AL (NK ) =XLT 
w a (NK ) =XK 
I F (S-APC) 16,41,41 

1(S ir (AHSF ( X-Xf^FAN) -Pnx ) 37 ♦524»b?4 
37 IF (AHSF (Y-Y'"'FAM)-PnY)40*524,52A 
4 0 rnriTiNiiF 
4] MHINT 17 

]7 format (//// 3nXl4HGE0nESTC ROUTE) 

HHTNT IR 

18 FOPMAT( 5X4hnAYS6XSHl.OMGI4X5H| ATI-5X4 HWAv/E 5X14H>*/AVE H I RECT I ON/2 X 9H 
1 C.F TRAVFL4X5H-TUnE/fX/^HTUnF5X6HHE:iGHT(SXioMFROM NORTH/) 
i: r?TNT ) q ♦ ( TF S ( K ) , A 7 ( K ) ♦ aL ( K ) ♦ H T ( K ) « W a ( K ) , K= 1 ♦ MK ) 

IH format ( F9,?,Fll,l,FR,l,F10,?,F13.0 ) 
r WE WILL NOW (TUFSS AN INITIAL VALUE OF THf CONTROL U(T) AND RUN A PATH, 
rjFLTAXrXFNO-XSTaPT 
i)FL. TAY = YFNn-YSTART 
.fF(l)F.LTAy .FO. o.O)GO TO 530 
aRGsOFI TaY/OFLTAX 
4| F ArAT AmF ( APG) 

]F(DEL7Ax)53l,530»53? 

53? IF (DELTay) 533»534,534 
633 v=6,2R31W ♦ ALFA 
00 TO 535 

534 v= A|.Fa 
GO TO 535 

531 v=3,l415R + ALFA 
GG TO b35 

530 IF (UELTAy) 536, 637,538 

536 V=4,7124 

(iO To 535 

538 v=l,57ofi 
GO TO 535 

537 PRINT 539 

538 F OrmAt (5xAHST 0P) 
f-0 TO P3 

535 CONTTMUF 
MMro 

GANA = f>.n 
541 LP=1 

r<X=FLOATF (KX) -I .01 
PKY = FL0ATF (K Y) -1 ,01 
M = 0 
L = 0 

0F| UMXxO,4 
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' D=.30 ,(I 

^ouDsionno . 
rftll=KTFS 
(>•■' boo 1 = 1 . ? 0 O 

^flll ll(I)=\/ 
bOl TI.AM(1)=1.0 

xM()( 1 )=n.n 

X'=^IG^’A (1 ) =(I.(I 
1I.AM(?)=0.0 
XMU(2) =1 .0 
XSTbWA (?) =0 . 0 
1 1. A M ( 3 ) = 0 . 0 
X'"(i ( ^ ) =n . 0 
XSIGMA (3) =1 . 0 
xsXSTAPT 
y=YST 6MT 
STF.P=TAU/?no, 

T=Hk/2A , 

f-'fl 50? T = 1,?00 

call tf.pp 

aPG=GAMA*M 
tall 3HTp 
K = I M I ) 

CAUL VDKPjv/ 
on 550 j=lt3 

XI AM ( J) = T|. AP ( J) - (CaPVX^COSF (II ( I ) ) *TLAM ( J) ♦CaPVX<»SINF (U ( I ) ) *XmII ( J) 
I ♦GAMA^^FvPF ( APG) *Hx<»x5TGMA ( J) ) *STFP 
xMi) ( j) =XMO ( j) - (capvy^cosf (u ( t ) ) «tlam ( j) ♦capvy«sinf (U ( T ) ) »XMIJ ( j) 

1 ♦GAMA*FXPF (aRGI-^MYAXSIGmA (J) )«STFP 
55 0 Tl AM ( J) = xLA^i ( J) 

XrX4.CAPV«>COSF (IJ ( T ) ) ^^CTFP 
Y = Y + CAPV<»S1NF (IM I ) ) <^STFP 
IP' (X ,L T, ?,nl ) GO TO 5?A 
IF ( Y ,LT, ?♦ 01 I PiO TO 5?A 
IC(X ,gT. PKX)(;0 TfT 5?4 
JF(y ,GT, PKy)GO to 524 
T=T+STFP 
502 COMTlNt'F 

Tmik 1 ) sxM|i( 1 ) 

T'‘ni (?) sXMiJ (2) 

flnjO('T= XLAM(n*XMi)(2)-XLAM(?)*XMll(l) 

U-AM( 1 ) =xMl)(?)/AD JOF7 
XMU ( 1 ) s-XMl) ( 1 ) /AnjDF.T 
x5IGma( 1 ) =0.0 
n AM (?) =_xLAM (?) /ADJDFT 
«f'ii (p) =X| AM ( 1 ) /aojoft 
X^^IGMA (2) =().u 

TLAM(3) = ( XI Am (2) «x MU (3) -XL AM ( 3 ) *TM|J ( ? ) ) /ADJDFT 

XMU (•<) = ( XL AM ( 3) *TMU( 1 ) - XLAM ( 1 ) *XM|J(3) ) /ADJDET 

XSIGMA (3) si .0 

XsxSTAPT 

rsYSTAFT 

7 = 0.0 

T=hR/?4, 

XTF5=0.0 

711=0,0 

712=0.0 

713=0,0 
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/?£> = (■) • 0 
Z?3=0.0 

733=0,0 

on 503 1 = 1 »«;oo 

CALL Tfpp 
M»G = nAVIA#H 
HT ( 1 ) =H 
TF<^( I)=XTFS 
CALL AMfil.F 
«7 ( I ) =xi n 
AL ( I ) =X| T 
w A ( I ) = X K 
call ShTP 
«; = I'd) 

CALL VDEPI'/ 

G»ADx ( I ) sI’LAM ( 1 ) * (CAPvd«COSF ( II ( I ) ) -CaPv*SIMF ( IK I ) ) ) ♦ 

1 (CAPVli*<;iK'F (IKI) ) ♦ CaPV/*COSF (0 ( T ) ) ) 

fWADY n ) =Ti AM (?) * (CAPVU«»cOSF (U ( I ) ) -CaPV*CINF (U ( I ) ) ) XMIK2)* 

1 (CAP\/|i*GTL'F (|)( T) ) ♦ CaP\/»COSF ( U ( I ) ) ) 

GRttU7 ( I ) =TLAM (3) * (CAPVlU^CnSF (U ( I) ) -CaPV«SINF (U ( I ) ) ) ♦ XMU(3)» 

I (CAPV/lioSlMF (IK T ) ) ♦ CaPV*C0SF (U ( T ) ) ) 

71 l=Zl 1 ♦n^AT^x ( I ) «GPADx ( I ) «STFP 
7) ? = ZlZ + (^PA['X ( I ) *0RAnY ( n «STFP 
713 = Z13 ^gPaDX ( I ) <»r?PAn7 ( I ) oSTpP 
722 = z22 + GP'M.)y ( I ) *GPADy ( I ) »STFP 
''?3 = Z? 3 *PRA 0 y ( I ) *GRAnz ( I ) »STFP 
733=733 *gPA07 ( I ) »gRAD7 ( I ) «STFP 
no 551 J=l,3 

xLaM( J) =TLaN’ ( J) -(CAPYX^COSF (II ( I ) ) <>TLaM( J) ♦CAPVX^STNF (U< I ) ) *XMU( J) 
1 4GAMA*FXPF ( ARG) ^^HX<^XSIGMA ( J) )*STFP 
XMli ( J) =XMil ( J) - (CARVYO'CnSF (U ( T) ) *TL.AM { J) *CAPVY*SINF (U ( I ) ) »XMH ( J) 

1 4GAmA*FXPF ( ARG) *HY*X5I6MA ( J) ) <>STFP 
5*^1 n AM ( J) =XLA-^ ( J) 

*=X*CAPV*COSF (U ( I ) ) «5TFP 
Y = y*CAPV»SIMF (IJ( I) )*5TFP 
Z=Z ♦ FXPF (ARG) *STEP 
1=T*STEP 
kTFS=XTF«;+GTFP 
503 COMTINUE 

xnoT = CAPV<»CO.SF (U (200 ) ) 

YOOT=CAPV«5ImF (U(200) ) 

ARG = GAPA<»H 
700T=EXPF ( aPG) 

CALL TEPP 
or (201 ) =M 
TFS (?01 ) =XTFS 
• CALL AN'GI F 
A7 (201 ) =xl.G 
Al (20l ) =XLT 
A, A ( 2 0 1 ) = X I' 
xo iFFsXF^'n-x 
Y0 IFF=YFmD-y 
GFLX=X-Fnx 
nrLYsY-FOY 

mSF= (DELX*<*24.oflY«*2-^973,75) / I 0A3,g 3G7A3 
FRR0P=xr)IFF«*2 + Y0IFF«»? 

IF ( FPPnR-M5F*MSF*CRIT»C9IT*0. 00002366 ) 505*505*511 
511 Al l=7ll*M>2*7l2*#?*xnnT<H»? 
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A1 ^=71 i<»71?*Z22»712*xDOT*YIJ0T 

tt2?»Zl2*«2*Z22«*?*YnOT**2 

f.'p’TsAl 1*A22-Al^««2 

C;i = ( A2?*xnTFF-Al 2»Yf)IFF) /OET 

C?=(All*Y0TFF-M2*xniFF) /OET 

e 1=Z11*C1*712*C? 

f- 2 = Z12#C1 ♦722*C? 

iiF|. T A1I = C 1 <^XI)OT*C2«YDOT 

1 AU = 1 AlunFl.T All 

'in boft K = l*2no 

uKl.U (K)=Fl*C-RAnX(K)+F2<^GPAF)Y(K) 

SOf, II (K ) =U (K ) ♦I’lELU (K ) 

|M = N* 1 

if(n .C’T. 2b)r,n to sqi 

i-iO TO 50 1 
so 7 PRINT 510 

bid F mpmaT (/bX 14HM0 CONVFt^GFNCE) 
bO TO rt3 
SOS comtinue: 

L=L> I 

I F (NM.E 1 ) GO TO S44 
IFHAi) .IT. TAyOuniOn TO S12 
OFl U^X=PE:U)^'X/?, 

IFlUFLI'i^X ,LT. 0,01) GO TO 514 
SI? TAUOLOsTAlJ 
fiO TO S4iS 

544 IF (Z.LT.ZOI D) GO TO S45 
l;FI UMX=OFU)MX/2. 

IF (DFLu^'X.LT.O.OI ) GO TO 514 

545 70l.0=Z 

Sas Kll= 7 ll* 7 ll ♦ 712*712 ♦ 7l3*7l3 ♦ XDoT*XOOT 

Hi2 = 7l?*7ll ♦ 712*72? ♦ 723*713 ♦ xi)OT*YOOT 

^13=2l3*Zll ♦ Z23*Zl2 ♦ l33*Zi3 ♦ XDoT*ZOOT 
h 22=712*712 ♦ 722*722 ♦ Z23*7?3 ♦ YnoT*YnOT 
823 = 713*71? ♦ Z23*Z22 ♦ 733*723 ♦ yIJoT*ZOOT 

G33 = Z13»713 ♦ 723*723 ♦ 733*7.33 ♦ 70oT*ZOOT 

i) = r1 1* (F)22*b33-h23*h 23) -Hl2*{Rl2*rt33.H?3*Rl3) 

1 ♦fU3*(H]2*B23-B22*Hl3) 

01 = (Hl3*ii?2-Bl2*R23) /D 
U?= (811*823-813*81 2) /n 
03- («l2#b 12-81 1*822) /L) 

H1=Z11«D1 + 712*02 ♦ 713*03 

S2=Z12*01 ♦ 722*02 ♦ 723*03 

83=713*01 ♦ z23*02 ♦ ?33*03 

b4=XD0T*0l ♦ YOOT*D2 ♦ Z00T*03 

8IG=0.0 

00 515 1=1,200 

OFLU ( I ) =P1*GRADX ( I ) ♦ R2*6RA0Y(I) ♦ R3*GRADZ ( I ) 
r)R = A8SE ( DEUKI) ) 

IF( BB ,GT, HiG )hIG=BB 
515 CONTINUF 

OR = OELUMX/( BIG ♦ 1 . o*ARSF ( R* ) /T aO ) 

00 520 1=1,200 
oFLU(I)=PR*DELU(I) 

520 0(I)sU(*I) ♦ OELU(I) 

IF (L ,6T. 25)60 TO 521 
GO TO 501 
5J4 CONTINUF 
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IF( MN ,GT. n ) fiO TO 54 0 
WPIM 5?? 

hf’P forma T ( ////2R X 1 flHM IN IMIJH TIMF ROUTE) 

PRINT 54.1 ♦» 

1^41 FORMAT ( //20X7HDET R sF6.?) 

RPINT Jfl 

PRINT IR. ( TFS (K) ♦ A7 (K ) »AL (K) ♦HT(K) ,WA (K) »Ksl »201 ♦ 10 ) 

OAmAsO , 0 1 
MMs 1 

r,r> TO 54 1 
540 PRINT 542 

542 format (////?«X iQHWFlOHTEn TImr ROUTE) 

PRINT 543,0 
PRINT IR 

PRINT 19» ( TFS (K ) ♦ A7 (K ) » AL (K ) tHT (K) ♦ WA (K ) »K = 1 *201 ♦ 1 0 > 

I’O TO P3 
5?1 PRINT 523 
523 format (5XPHMO SOL N) 

00 To Hi 

S?4 PRINT 525. <*Y 

^;,5 f ormat (5X4f|HSHIP TS MOW OUTSTOE ThFATEP OF OPFPATION 5 X 3 HX =Flo.5» 
1 5x3hY spin, 5) 

HO STOP 
F . \in 
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subroutine TERP 

DIMENSION C(4)»P(4)fO(4),PX(4)*QY(4), HP (4)» CP (4)fSP (4), HP X(4), 
1hPY(4) ,CPX (4) fCPY(4) ♦SPX(4) fSPY(4) ,HO(4) ,00(4) ♦S0(4) ,HS(4) ,CS(4) , 
2SS(4) ♦HT (4,4) ,CT(4,4) ,ST (4,4) 
common X, Y, AtR*CC*H*CK«SKtKX*KY,Kl,FOX,FOY,KcOM 
COMMON/L l/XHT (fl44fl) ♦CSK (fl448) ,SNK (8448) »T 
COMM0 N/L?/HX,Hy 
COMMON/L l/DKX *DKY 
L = XF TXF(T) 

TT=(-INTF(T)^T)«2.-1. 

TP1=IT+1 . 

T^n=TT-1 . 

T?m=TP1«tmi 
IF (L-KI^3) 1.1,4 
I IP(L)?,?.3 
? K4 = 3 

Tm3=TT-3. 

CM) STM10TM3/H. 

C(2)=-IP1«tn3/4, 

C (3) =T?M/fl. 

00 TO 16 

3 K4s4 

G= ( J.*TT4>2. ) *TT-R. 

F = -4,*TI.r, 

C ( 1 ) =-T?m«Tm]/i6, 

C(2)=G<»TM1/16, 
r (3 ) =-f»tp1/16. 

C (4 ) aT2M»TPl/16. 

GO TO IB 

4 IF (L-KT)6.5*5 

5 k4s1 
LsKI-1 
C(1 )=1. 

GO To l6 

6 Gs(3,*TT*2.)*TT-9. 

K=-4 .^TT+G 

C (1 ) =-T?m*Tm1 /16, 

TF(L-KT+1) 11,10,5 
in k4=2 

C (2) = (G*TMU (T2M.F) *TP1 ) / ) 6, 

GO TO 15 

11 C (2) =G*TM1 /16, 

IF (L-KT^?) 3, 12,10 

12 K4=3 

C(3) = (T2^^-F)«TP1/16. 

15 L=L-l 

16 m=XFTXF(X)-2 
N=XF IXF ( Y)-2 
x,y=(-imtf(X) ♦X)»2,-l. 

YY=(-IMTF(Y)+Y)#2,-1. 

XP1=XX+1 , 

XMlaxX-1 , 

YP]=YY^1, 

YMi=YY-1 . 

X2msXP1«XM1 

Y?MsyP1*y^^1 
P ( 1 ) =-XMl *X2M/16, 
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P(?)=((3,*XX*2,)*XX-9,)*XMl/lft. 
P (3) s-XX*XX/8, +1 , 125-P (2) 
P(4)=XP]*X2M/16. 

0(1) =-YPl »Y2M/lft, 

0(2) = ((3,<»YY*2,) h»YY-< 9. ) *YM1/16. 
iJ (3) s-YY<»YY/fl. + I .125-0(2) 
ij (^ ) =YPl *Y2 m/1^, 

)7 PX (4)e(0,375*XX-0.125)*XPl 
PX (1 )sO,S*XX-PX (4) 

PX(2)s(1 ,l?!5*XX-1.375)<tXpl 

PX (3) =-n,5«XX-PX ( 2 ) 

OY (4)s(0.375*YY-0.1?5)*YPl 
uY { 1 ) =0,5#YY-(JY (4) 

Oy (2) = ( 1 , 125»yy-1 ,375) <>YP1 
OY ( 3) =-0,5*YY-OY (?) 

Itt DO 27 Ksl *K4 
HP ( K ) *0 , 0 

ro (K) so.o 

SP(K)s0,0 
10 HPx(K)s0.n 
HPY (K) =0.0 

cpx (K) =o,n 

('OY (t'')=0,0 
sox (K-) =0,0 
SPy (K) =0,0 

20 K>^=( (K>D<»KY^N )«KX 4. M -KC(0M 
C)0 23 J = 1.4 

HO ( J) -() , 0 
CO ( J) =0.i> 

SO ( J) =0.0 

21 HS(J)=0,0 
CS ( J) =0 .0 
s*; ( J) =0.0 

?? JJ=J*KX*KK 
23 1=1 t4 
[ I=I ♦ JJ 

•^Td.jisxprdi) 

CT ( I . J) =CSK (IT) 

23 ST ( 1 , J) =SMK ( TI) 
iio 25 1 = 1.4 

no 25 J=] .4 

nn ( I ) =0 ( J) WHT ( I . J) ♦Hf)( I ) 

Cl) ( I ) =0 ( J) *CT (1 , J) >CD( T ) 

50(1) =(J ( J) «ST ( I . J) ♦SO ( I ) 
HS(I)sP(j)»hT(j.d+HS(l) 

CS ( I ) =P ( J) *CT ( J. T ) ♦CSd) 
S<?(I)=P(J)<»ST(J.I)>SS(T) 

25 cowl I Nil ip 

00 27 1 = 1 .4 

HP (K) =HO (I ) *P ( I ) +HP (K) 

CP(K) =CO ( I ) ( I ) +CP (K) 

SP(K)=SO(I)*P{I)+SP(K) 

?H hdx (K) =H0 ( I ) *PX ( I ) ♦HPX (K) 

CPX (K) =cn ( I ) *PX ( T) +CPX (K) 

SPX (K) =S() (I ) *PX ( I) +SPX (K) 
mPY (K ) =HS ( I ) «QY ( I ) ♦HPY (K) 
CPY(K)=CSd)*QY(I) *CPY(K) 

SPY ( K ) =SS ( T ) *(1Y ( I ) +SPY (K ) 
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?J COIMTINIIF 
HsO • 0 
CKaO.O 
SK=0.0 
HXrO.n 
HY = 0 , 0 
CKX=0.0 
CKY = 0. 0 
•^KXsO .0 
'^KysO.O 
(10 31 K = 1,K4 
H=C (K) «HP (K) ♦H 
CK=C(K)«CP(K)*CK 
SK=C (K) «SP (K) 4.SK 
OX=C(K) <»mPX (K) *HX 
Hy=C(K)*hPy(K) ♦Hy 
CKX = C(K)<»CPX(K)*CKX 
C<Y = C (K) <»CPY (K) *CKY 
SKX=C(K)*SPX (K) ♦SKX 
Sky = C(K)<»SPY(K)*SKY 

■^1 COMTINUf: 

KADrSOPTF ( CK*CK*SK*SK) 

rK=CK/PAD 

SK=SK/RAn 

32 rjKX = CK«SKX-SK*CKX 
r)KY = CK*SKY-SK*CKY 

3’’ WFTUPN 



SI'F-ikOUTlMf- 

common Y.o , 8,00* H,CK«SK»KX» KY,KTf fox tFOYtKCON 

rOMMON/L ♦ HY 
COMMON/L4/aX, AY»BX,RY»CX,CY 
SPFP.n PAPAMF^TFOS of a API 
^ n OBHRn 7*^ 

vH = FACl«*?. I A46Q/fhBRH«FXPF ( -0 . OSi 3470S93044*H ) *0.0253410549^3 
(■ '^Clr2, 1 49i494f,?5g/pACl 
F AC2=FAr) -0.093470593044 
nVH = VH*F aCF^ 

^4C1=H*H.045AH6?274 

vi- =F AC l**(i. 34H3 12 j959B«typF (-n.025324 954736*H) *7.9637908524 
F AC 1=0. 34 H n?l9696/FACl 
F AC2 = FACl-() .0 25324954736 
In'F =\/F*FAC2 
h ACl=H*37. 097335211 

M = F An **3,58243605 31 *FxPF ( -0 . 0 88 02945249 0 *H ) *0 . 40 R48F?7329 1 F-4 

F AC 1=3.5 w 24 36 058] /FAC I 

t-AC2 = FA r 1-0. 088 0 294 52 490 

3(< = ri*F A02 

Q= ( VF ♦ VH) *ii .5 

CC=A-VH 

i »A= ( O'v/F ♦OVM ) *0 . 5 

l)C = UA-f)VM 

A X = Ua*F^X 

aY=UA*Hy 

hX=Q8*mX 

‘'Y = D3*HY 

C Y=DC*hV 

CY=l)C*hy 

‘-'c’TUPm 

F MO 
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SHRKOUTINE VOERIV 

COMMON X f Y* A *0 tCC »H tCK ♦ SK t K X « K Y t K T ♦ FOX » FOY » KTON 

C0MM0N/L3/l)KX»DKY 

COMMON/L4/AX , AYf BX»BY»CX tCY 

C0MM0N/L5/CAPV»CAPVX»CAPVYfCAPVUf W 

REAL MSF,MSFX»MSFY,NUM 

DF.LXaX-FOX 

ofly=y-foy 

MSF= (OEL x*»2 + DELY**2 + 973.75) /1043.63P743 

MSFx»nELx/521 .0193715 

mSFY=DELy/521 .0193715 

IF( LP •kQ* 0 )60 TO 1 

C.OST=C05F{W) *CK ♦ SINF(w)*SK 

STNT=SIMF (W) *CK - COSF(W)*SK 

GO TO ? 

1 COSTsCOSP*CK ♦ SINP*SK 
STMTsSINP^CK - COSP4SK 

2 a2mC?=A*A-CC*CC 

POOTsSQRTF (B*«2* (COST»*2) * A2MC2* { 0 I NT**2 ) ) 

Ml im = R*A2»iC2/a 

UFNOM=B*CC^»CnST/A*ROOT 

v=num/DEmom 

VGWiOaV/H, 56604 10667 

CAPVsVGRin*MSF 

amCx=A*Ax-CC*CX 

AMCY*A*Ax/CC#Cy 

FX= (CC*PX+0*CX-R*CC*AX/A) /A 

FY= (CC«PY^R*CY-R*CC<»AY/A) / A 

BPCMA2sR«R4.CC#CC-A*A 

TERMX* (C0ST*FX^84CC*SINT*DKX/a+ (AMCX*SINT*SINT+b*RX*cOST*C05T ♦ 
1 RPCMA2*SlNT*c0ST*DKX)/R00T)*(A*V6RIn/(B*A2Mc2> ) 

TERMys (C0ST*EY4.B*CC*5INT*DKy/A+ ( AMCY*SINT*SlNT + R#PY*rOST*COST 
1 ♦8PCMA2*SINT*COST*DKY)/ROOT)*(A*VORID/(R*A2MC2) » 
vxsVGRID* (2.*amcX/a2MC2 ♦BX/0-aX/a-TfRMX) 

V Y = VGRI^'^^ (2,*AMCy/A2MC2 + By/B-Ay/A-TEPMy) 

CAPVXsVX«MSF+VGRID*MSFX 

CAPVYaVY«M5F^VGRlO*MSFY 

VM-VGRIO* (B«CC*SINT/A + 8PCMA2#SINT*CO‘;T/ROOT) / (B#CC*CoST/a +ROOT 
r,APVLIsMSF*VU 

return 

F ND 
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■I 



c;iiPHOUHNF flMGLF 

COMMOM X«Y«A»H»CC»H«CKf SK » K X ♦ K Y ♦ K T ♦ FOX » FO Y t KCON 
COMMnN/t.6/ XKiXLGtXLT 

I X = X-FOX 
,',FL Y = Y-FOY 

rojxKs-oi- L x«rK-nFLY»c;K 
KTMXKsOFI,X»«;K-[)FLY*Ct< 

IF (COSXk ) ?» 1 «2 
1 «K = bTF;MF(qO.»SINvK) 
sn -n (S 

? xKsATANF (9IMXK/Cn<;XK)»57. 29577951 
TF (COSXK) 

3 IF (STMXK ) 5,4,4 

4 xk = xK4.iho, 

00 TO 5 

5 xK=x'^-i^n. 

^ r F ( X K ) 7 , F' , « 

/ XKrJOO.+XK 

xT=UFLx<>,9<i4Hn775-nF.|,.Y<>. 17304518 
VT = UFI x<f,i 73048] «*DELY<*. 98480775 
waO=SOKTF ( xT«XT ♦ YT^t’YT) 

1 F ( X T ) ] 0 , 9 , 1 0 

') XI, g = SIONF (90. ,YT) 

■Hn |o ) <*• 

l!) x|.6 = ATAN'F (YT/XT) o5 7, 2957 795] 

T F ( X T ) n , 1 4 , 1 4 

11 IF(YT)13,12,)? 

12 xL(-. = XLG 15(1. 

10 14 

] 3 xLG=xL(’ - 180. 

]4 XI T=-ATAnF (RAO/ 31 ,205) *114.591559 ♦ 9o, 

RfTURM 

FM|' 
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TDENT DTGB 

PPOGRAM LFWGTH 
BLOCKS 

t)ROGPAM<> LOCAL 



ENTRY POINTS 



oonouo nTGfl 
F.XTFpmAL symbols 
TBPOTG 



« 

* 

« 



PTGn 



entry otgb 

FORTRAN CALL 
---ydHERE A ■ 
B B 



TO THIS program-— CALL DTGB(A,B,C) 
address of DTQ W0RD(LEFT JUSTIFTEn)Bl 
ADDRESS TAu B2 

address new DTG W0PD(LEFT JUSTIFIF0)B3 



return 


with mew dTg in C 


BSs7 


1 


sxft 


H3 


SAiS 


T 


SAi 


Bl 


l-Xl 


A a 


SA? 


R2 


rj 


sXTBFDTG 


SA 1 


T 


SB^ 


XI 


LXG 


12 


saa 


R3 


JP 


DTGB 


data 

end 


0 



RGT adjust DTG 



UNUSED STORAGE 



21 statements 



3 symbols 
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TDPNT FIELPSX 
5P0GPnM lfnGTH 

51.0CKS 

LnCAL 



"MTPY POTmTS 

oonnoo fifldsx 



■XTERNAL SYMHOLS 



ryOn T wAp64X RDmST 



0 

» 

» 

» 

If 

If 

» 

If 

If 

» 

If 

If 

F IFl. HGX 

If 



If 

BGM 



» 



If 



F^'TPY FIELDSX 

oFvISEn BY SEl fridge JUM IRGG MELLONICG 

GEl.ECT FNWF FIFl.D FPOM MASTER* UNPACK AND FLOAT DATA 

FORTRAN CALL TO THIS PROGRAM-- CALL F lELDS ( A » P t C » D) 

WHERE A B ADDRESS OF 39BR WORD RINCCONTAlNS FLTD DATA ON EXIT) 

fl = address of date-time GROUP WORD 
c = address Of Parameter selection index (o»i #2*3) 
where 0 s CD#l s CH»2 s u CURR»3 * V CURR 
D S address of read-master error indicator word 

<tf fteLo not FOUNp OR PaRITy error (d>*i on Return) 



HSS 7 1 

SAv/E Fortran call parameters 

SA? 92 

SX6 rtl 

RX7 X2 

<;Ais param 

SA7 TABLE+1 

SA3 R3 

HXH X3 

SA(S COUNTER 

SX7 B4 

SA7 PARAM+1 

pead One Field from master tape 

SA2 counter 

SAi NAMTAB+x? 

HXf, XI 

SAft TABLF 

SHI A6 

RJ xXROMST 

unpack field 

G AO rO 

SHl RAF 

SA? PARAM 

SB? X2 

SB3 7601B 

RJ *XWAB64X 

shift data into lower Afi BITS 

SBl 20 

GB? 3908 



OTG 

STORE selection INDEX 

not-found indicator 

BCD parameter name OF FIELD 

Loc. unpacked field 
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SB3 


BO 




SB? 


1 




SAl 


PARAM 


SHIFT 


SA? 


XU01 




AX2 


12 




RX6 


X2 




SA6 


Xl*03 




SB 1 


01*B7 




SB3 


B3*B7 




GE 


B2»B1#SHIFT 


* 


convert field to floati 




SBl 


Xl 




SB2 


XI 




SAl 


counter 




SA2 


SCAtAB*x1 




SB3 


X2 




SB4 


7601R 




RJ 


■XFXDFLT 




JP 


fieldsx 


TAPI F 


data 


0 




data 


0 




VFO 


30/PAF,3o/ErR 




data 


5LTAPE4 




DATA 


3LEEE 


NAMTAP 


data 


lOHCD 




DATA 


lOHCH 




data 


lOHU CURR 




data 


lOBV CURR 


SCAT ah 


data 


41 




data 


41 




data 


36 




data 


36 


COUNiTFR 


[)AtA 


0 


PAPAM 


BSSZ 


2 


ERR 


SAl 


PARAM*1 




SBl 


1 




SX6 


01 




SAft 


XI 




JP 


FIELDSX 


PAF 


BSSZ 

ENO 


1100 



UNUSPn storage 85 STATEMEmTS 



scale factor 



13 SYMBOLS 



THTMT SEALAMD 
PPOG»AM l.pK'GTH 

^-LOCKS 

POOGPAM# LOCAL 

FMTPy points 

OflOOOfi c;EALANn 
txrspMAL. «;YMOnLS 
lani)«;f A 



cNtRY SEALAMD 





FOoTPAN CALI. TO THIS PROGRA 


M — - CALL SEALANO ( I » J»LaND»KER) 


o 


-WMfPE 


I * ADDRfSS FO RfaL 


T index 


ON FNWc 63X^s3 GRId 


o 




J s ADORFSS of RFAL 


J index 








LAMO = address of LaND/SFA 


indicator (0 FOR SEA# 1 FOR 






KER = address of error indicator (OUT OF BOUNDS) 






(0 IF ERROR absent* 


1 IF error is PRESENT) 


A Ml) 


MSS7 


1 








SA] 


R1 




I 




SA? 


P2 




J 




HXl 


XI *R1 








IX? 


x?*r2 








SHl 


R1*3R 








SH? 


B2*15 








1. XI 


X1*B1 




1 S39 




l-x? 


X2«R? 




J Sl5 




MXO 


36 








.4X1 


xo*xi 








HX? 


-X0»y? 








Bxi 


xi*x? 




1 BITS 24.47»J BITS 0-23 




SXft 


Pi 




EACH SCALED Sl5 




SX7 


P4 








SAft 


LSIM 




SAVE L/5 TND. ADDRESS 




sat 


FRPIN 




Save error ind. address 




SB], 


SEA 








SHp 


land 








SH3 


error 








,|P 


sXLANDSFA 






SLA 


MX A 


0 








MX7 


0 








,|P 


OUT 






LAMP 


sxa 


1 








MX 7 


0 








JP 


OUT 






FWPor^ 


SXa 


1 








SX7 


I 






OUT 


SAl 


LSIN 








SA? 


ERRIN 








SA6 


XI 








SA7 


X2 








jP 


SEA|_amD 







LANf, 



LSIN 

ERRIM 



HSSZ 

KSSZ 

FMO 



JMiMr IP 



OIMFMSIoN of wavF FIFLD AMRayS=( 32t 12) T-iU'= In 

WOUTF OF '^HIP HF.GINS 0 HOURS AFTPP TIME SERIES 0«TGIN 
FROM UOMOITHDE = 154. n ANM> LATITUDE = 41.0 

TO LOMOITUOE s -i?3.n <\Nn LATITUDE = 3P.0 



OA YC 


long I 


TRA VFI 


-THOE 


n , On 


154.0 


.5r, 


1S5.5 


1 . On 


15Q.9 


1 .5f) 


lH3.t3 


2 , On 


1H7.S 


2.5r 


171.0 


3.00 


175.3 


3,5o 


179,6 


4.00 


-1 /5.4 


4 .So 


-170.3 


5.00 


-165.2 


5.5o 


-160,0 


H.OO 


-155,0 


H . 5^ 


-150,1 


7.00 


-1A5.S 


7. So 


-141,2 


5 .00 


-135,8 


5. So 


-132,6 


R.no 


-125,5 


R.sn 


-124,6 


9,7n 


-123,1 



OFODESTC 



1 ati- 


wave 


TUnt, 


HEIGHT 


4 1,0 


?l ,21 


4 2,1) 


19,31 


4 3 , 0 


16.76 


44,2 


15.63 


45 , 1 


20.00 


45 , H 


19.54 


45,5 


16.96 


47.1 


1 1 ,64 


4 7,5 


9.71 


47.7 


10.36 


47.7 


4,09 


<»7.4 


6.55 


A 7 . 0 


9,48 


46,3 


11.78 


a5.4 


17.28 


44,4 


13.41 


43.2 


7,39 


41.6 


4,86 


40.2 


1.91 


38.5 


6.44 


37.8 


13.70 



Wave niREcTiON 
from north 

3l R 

3n4 

231 

25? 

27 f) 

IRR 

2F-5 

?1R 

lOR 

20 H 

121 

23? 

150 

175 

175 

15 ? 

153 

162 

P9 

345 

314 



45 



MIMIMIIM TIMF ROliTF 





DEI 


M * 


,00 




oayr 


LOMOI 


1 ATl- 


mAvE 


WA\/E DIPEcTIOn 


OF TPavfI. 


-Tlt()E 


TUnE 


height 


from •'IOoTh 


0 . 00 


154,0 


4 1.1' 


21.21 


319 


UP 


157,7 


4 0,4 


14.80 


1 


.94 


IM .b 


4 0.3 


6.92 


171 


1 .44 


166,2 


40.5 


10.10 


291 


1.9? 


170,2 


40.9 


14.70 


159 


2.40 


174.3 


4 1.6 


15.68 


?8l 


2. ftp 


178.1 


42.7 


17.50 


28 0 


3.36 


-17H.0 


44.0 


13.88 


275 


3.84 


-173.8 


45.1 


12.69 


282 


4.3? 


-169,3 


45.9 


12.48 


275 


4 ,Ro 


-164.7 


46.0 


7.06 


142 


b.?P 


-160.0 


45.6 


6.89 


334 


5.74 


-155.6 


44.9 


9.35 


147 


6.24 


-151.3 


44,3 


10.78 


157 


6.7? 


-147.3 


43.7 


13.81 


164 


7.?n 


-14 1.4 


43.1 


12.97 


171 


7.6« 


-119.3 


42.5 


9.42 


166 


8.16 


-135.1 


4 1.7 


5.62 


296 


H .64 


- 130.8 


40.7 


3.96 


115 


9.1? 


- 1 26 . B 


39.4 


1.29 


1 


9,60 


-123.0 


38 , 1 


15.20 


8 



.,FIGHTEn TIME ROllTF 



OF 



pf:t » = .^2 



OAYc 


LONG I 


TRA\/p| 


-Hint 


O.On 


154,0 


.48 


157.6 


.97 


161,6 


1 .4? 


166,0 


1.94 


170.2 


2.4? 


174.2 


2.9n 


177,9 


3.30 


-178.4 


3.87 


-175,3 


4.36 


-169,7 


4.84 


-164,9 


5.3? 


-160.1 


5.81 


-155.6 


6.29 
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A, 



GLOSSARY 



TEST 



IMIN JMIN 




Minimum indices of wave heiqht 
data subfield 


KX, KY, KT 


= 


Dimensions of wave data sub- 
field 


DTG(KT) 


= 


Array of octal-date time groups 


U(200) , DELU(200) 


= 


Control arrays 


GRADX, GRADY, GRADZ 


= 


Gradient arrays 


AZ, AL, HT, WA 


= 


Arrays for track position 
tabulation 


WE(3989) , WD(63,63) 


= 


Unpacking area of core 


XLAM, XMU, XSIGMA 


= 


Components of adjoint vectors 

12 3 

, A , a"* 


X,Y 


= 


Rectangular coordinates of ship's 
position relative to subfield 
origin 


A,B,CC 


= 


Parameters for elliptical polar 
velocity diagram Fig. 1 


H, CK, SK 


— 


Wave height, cosine and sine 
of wave direction relative to 
OXY axe s 


FOX, FOY 


= 


Floating point coordinates of 
North Pole relative to subfield 
origin 


XHT, CSK, SNK 


= 


Wave field arrays 


T 


= 


Time in days from first member 
of time series 


TFS 


= 


Time in days from beginning of 
track 


xmean, YMEAN 


= 


Coordinates of midpoint of sub- 
field 


ROX, ROY 


= 


Semi-dimensions of a rectangle 



in the subfield within which 
track must lie 
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LR 


= Index to denote a geodesic track 
(LR=0) or an ootimum track 
(LR=1) 


COSLGl, SINLGl 


= Cosine and sine of longitude of 
initial point 


COSLTl, SINLTl 


= Cosine and sine of latitude of 
initial point 


COSLG2, SINLG2 


= Cosine and sine of longitude of 
terminal point 


COSLT2, SINLT2 


= Cosine and sine of latitude o'F 
terminal point 


XI, Yl, X2, Y2 


= Coordinates of initial and ter- 
minal points of track relative 
to North Pole in grid units 


S12 


= Straight line distance from XI, 
Yl, to X2 , Y2 


ARC 


= Great circle distance from XI, 
Yl to X2, Y2 


XSTART, YSTART 


= Coordinates of initial point in 
grid units relative to s\±ifield 
origin 


XEND, YEND 


= Coordinates of terminal point 
in grid units relative to sub- 
field origin 


COSP, SINP 


= Cosine and sine of ship's head- 
ing at any point along great 
circle track 


V 


= Control program to define 
initial control u(t) 


GAMA 


= Weighting factor used in cal- 
culating optimum track, GAMA=0 
minimum time, GAMA > 0 weighted 
time track 


DELUMX 


= Maximum allowable change in the 
control u(t) at any instant of 
time 


N, L 


= Counters to stop computation if 
convergence is questionable 
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STEP 


= Variable time steo of integra- 
tion 


XDOT, YDOT, ZDOT 


= Time derivative of state 
variables 


MSF 


= Map scale factor for stereo- 
graphic projection 


All, A12, A22 


= Elements of A matrix of Sec- 
tion III 


Zll, Z12, Z13, Z22, 
Z23, Z33 


= Elements of Z matrix of Sec- 
tion III 


DET 


= Determinant of A matrix of 
Section III 


Cl, C2 


= Elements of C matrix of Sec- 
tion III 


El, E2, DELTAU 


= Elements of E matrix of Sec- 
tion III 


NN 


= Index to denote a minimxim time 
track (NN=0) or a weighted 
time track (NN=1) 


Bll, B12, B13, B22, 
B23, B33 


= Elements of B matrix of Sec- 
tion III 


D 


= Determinant of B matrix 


Dl, D2, D3 


= Elements of D matrix of Sec- 
tion III 


Bl, B2, B3 


= Elements of a matrix of Sec- 
tion III 


RR 


= Coefficient of DELU(T) to limit 
maximum change in control at 
any time to DELUMX 

TERP 


C(4) 


= Weighting factors for interpola 
tion in the time demension be- 
tween K4 ordinates 


L 


= Index to determine member of 
field time series 


TT 


= T of Eg. (20) of reference [5] 
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M, N 

XX, YY 

P(l) to P(4) 

Q(l) to Q(4) 

PX, QY 

HX, HY, CKX, CKY, 
SKX, SKY 

DKX, DKY 
VH, VF 

AX, BX, CX, AY, BY, 

MSFX, MSFY 
COST, SINT 
VGRID 
CAPV 

CAPVX, CAPVY, CAPVU 



= Indices to pick out x, y grid 
point data 

= XJf of equation (18) of 
reference [4] 

= Pi to P4 of equation (19) of 
reference [4] 

T 

= Elements of P (y) matrix of 
equation (18) of reference [4] 

= Partial derivatives of P(l) to 
P(4) and 0(1) to 0(4) 

= Partial derivative of H, cos K, 
sin K 

= Partial derivatives of wave 
direction angle K 



SHIP 



= Ship's speed in knots in head 
waves and following waves 

CY = Partial derivatives of param- 

eters of elliptical oolar 
velocity diagram 

VDERIV 

= Partial derivatives of map 
scale factor 

= Cosine and sine of the angle 
0 of Figure 1 

= Ship's speed in knots relative 
to the stereographic grid 

= Ship's speed in knots relative 
to the earth's surface 

= Partial derivatives of CAPV 
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